
! This is for extracting maximum temperatures and fields from
! the full simulation runs for the purpose of photodisruption
! plotting in pd_plot.m
! SWH 2011

program extract_max

implicit none

real(8)       :: time, evar, pvar, wvar, efvar, envar
real(8)       :: emax, pmax, wmax, efmax, enmax
integer       :: temp_file=10
integer       :: rtemp=11
integer       :: tstep, ios, fnum, ii
character     :: junk
character(4)  :: cii

fnum = 4000

open(rtemp,FILE='../results/T_max.txt',STATUS='UNKNOWN')

do ii = 1, fnum

  if (ii .ge. 1000) then
     write(cii,'(I4.4)') ii
  elseif (ii .ge. 100) then
     write(cii,'(I3.3)') ii
  elseif (ii .ge. 10) then
     write(cii,'(I2.2)') ii
  else
     write(cii,'(I1.1)') ii
  endif

  open(temp_file,FILE='../results/t'//trim(adjustl(cii))//'.txt',STATUS='OLD')
  read(temp_file,*) junk

  emax = 0.0D0
  pmax = 0.0D0
  wmax = 0.0D0
  efmax = 0.0D0
  enmax = 0.0D0

  do
    read(temp_file,*,iostat=ios) tstep, time, evar, pvar, wvar, efvar, envar
    if (is_iostat_end(ios)) then
      close(temp_file)
      write(*,'(a,I5,a)') " >>>", ii," complete"
      exit 
    endif

    if (evar .gt. emax) emax = evar
    if (pvar .gt. pmax) pmax = pvar
    if (wvar .gt. wmax) wmax = wvar
    if (efvar .gt. efmax) efmax = efvar
    if (envar .gt. enmax) enmax = envar
  enddo

  write(rtemp,'(E20.12,1X,E20.12,1X,E20.12,1X,E20.12,1X,E20.12)') emax, pmax, wmax, efmax, enmax
enddo

close(rtemp)

end program extract_max
